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We obtain the complete phase diagram of the hardcore Bose-Hubbard model in the presence 
of a period-two superlattice in two and three dimensions. First we acquire the phase boundaries 
between the superfluid phase and the 'trivial' insulating phases of the model (the completely-empty 
and completely-filled lattices) analytically. Next, the boundary between the superfluid phase and 
the half-filled Mott-insulating phase is obtained numerically, using the stochastic series expansion 
(SSE) algorithm followed by finite-size scaling. We also compare our numerical results against the 
predictions of several approximation schemes, including two mean-field approaches and a fourth- 
order strong-coupling expansion (SCE), where we show that the latter method in particular is 
successful in producing an accurate picture of the phase diagram. Finally, we examine the extent 
to which several approximation schemes, such as the random phase approximation and the strong- 
coupling expansion, give an accurate description of the momentum distribution of the bosons inside 
the insulating phases. 
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I. INTRODUCTION 

One of the most remarkable achievements in the field 
of ultracold Bose gases in recent years has been the ob- 
servation of a superfluid to Mott-insulator transition in 
optical lattices i By playing with the intensity of the dif- 
ferent laser beams involved in the setup, experimentalists 
have been able to study this transition in effective one^ 
two^ and threei dimensional geometries. This extraordi- 
nary accomplishment was achieved with gases of bosonic 
atoms confined in optical and magnetic traps. Using the 
strength of the optical lattice MS c\ control parameter, 
these gases were reversibly tuned from a Bose-Einstein 
condensate to a Mott insulator (a state composed of lo- 
calized atoms) X 

It is generally accepted that this quantum phase tran- 
sition can be studied using the Bose-Hubbard model, 
where the transition is found to be from a compressible 
superfluid phase to an incompressible Mott-insulating 
one (SF-MI)£ Over the years, much theoretical work 
has been devoted to determining the phase diagram of 
the model in various dimensions, using many different 
approaches However, a direct comparison between 
theoretical results and experimental ones^— still remains 
obscured by issues such as the spatial inhoniogeneityj 1 ^— 
finite-temperature effects ) 16 ' 17 and the limited set of ex- 
perimental tools available to probe the nearly isolated 
ultracold atomic systems. 

In a recent paper, Aizenman et alr^ argued that 
the phases of the Bose-Hubbard model can be studied 
equally-well by examining a slightly different variant of 
it, namely the Bose-Hubbard model in the limit of infi- 
nite onsite repulsion (i.e., the case of hardcore bosons), 
in the presence of an alternating (checkerboard) onsite 
chemical potential (a superlattice with period two) . The 
advantage of studying the latter model lies in the fact 



that it exhibits all the salient properties of the Bose- 
Hubbard model, while also being more amenable to an- 
alytical treatment. Specifically, Aizenman et al. rigor- 
ously proved the existence of SF and MI phases in the 
half-filled three-dimensional case (although they did not 
show that there is no intermediate phase between the 
two). In Ref. [H two of us (I.H. and M.R.) studied that 
very same model for the case of zero chemical poten- 
tial both in two and three dimensions, using quantum 
Monte Carlo simulations and analytical approximation 
approaches. We showed that the SF-MI phase transition 
is a direct transition, and we determined its critical value. 

The hardcore Bose-Hubbard model with a superlat- 
tice has yet another attractive feature that the general 
Bose-Hubbard model lacks: it is exactly solvable in one 
dimension. This is due to the existence of a mapping 
of the hardcore bosons to noninteracting fermions. This 
in turn enables the evaluation of correlation functions of 
interest by exact means 

In this paper, we study the complete phase diagram of 
hardcore bosons in the presence of a superlattice in two 
and three dimensions and with arbitrary chemical po- 
tential. We determine the phase boundaries separating 
the compressible SF phase of the model from the various 
insulating phases. First we acquire the phase bound- 
aries between the SF phase and the 'trivial' insulating 
phases (the completely-empty and completely-filled lat- 
tices) analytically. Then we perform high-precision nu- 
merical simulations using the stochastic series expansion 
(SSE) algorith m 23 : 24 in order to find the phase bound- 
ary of the transition between the SF and the half-filled 
MI. This is done by calculating the free energy the 
density of bosons in the zero-momentum mode po, and 
the superfluid density p s . The latter two quantities drop 
to zero upon entering the insulating regime from the SF 
phase. 
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Once the complete phase diagram is obtained, we pro- 
ceed to examine the model analytically by employing two 
mcan-hcld-typc approximations and a strong-coupling 
perturbation scheme (up to fourth order in the hopping 
parameter) in order to determine the extent to which an- 
alytical methods allow a reliable description of the system 
and its various physical properties, specifically in the con- 
text of the phase boundaries separating the compressible 
SF regime from the incompressible insulating regions. 

The paper is organized as follows. In Sec. [TT] we re- 
view the model at hand and present a qualitative de- 
scription of its expected phase diagram. In Sec. IIII1 we 
compute the phase boundaries between the SF phase and 
the empty and filled lattices analytically. In Sec. |lVl we 
obtain the remaining boundary between the SF and the 
half-hllcd MI phase. This phase boundary is computed 
numerically, using the stochastic series expansion (SSE) 
algorithm. Section [V] is devoted to studying the phase 
diagram as it is given by two mean-held approaches, and 
in Sec. I VII we employ a strong-coupling expansion (SCE) 
method. These approximation methods are then com- 
pared against the previously obtained numerically-exact 
results. In Sec. I VIII we study the momentum distribu- 
tion of the bosons, in order to allow for a comparison 
with future experimental data. Finally, in Sec. IVIII1 we 
conclude with a discussion and summary of our results. 

II. MODEL 

The Hamiltonian for hardcore bosons in a period-two 
hypercubic superlattice in d-dimensions, with N = L d 
sites and periodic boundary conditions, can be written 
as: 

(ij) i 1 

(1) 

Here, (ij) denotes nearest neighbors, a, (aj) destroys 
(creates) a hardcore boson on site i, hi = djdj is the 
local density operator, fj, is the global chemical poten- 
tial, and A(— l) cr ( l ) is a checkerboard local potential with 
a(i) = on the even sublattice and 1 on the odd sub- 
lattice. The hopping parameter t (which we shall fix at 
t = 1) sets the energy scale, and without loss of generality 
we choose A > 0. The hardcore boson creation and an- 
nihilation operators satisfy the constraints dj 2 = af = 
and {di,aj} = 1, which prohibit double or higher occu- 
pancy of lattice sites, as dictated by the U — > oo limit 
of the Bosc-Hubbard model. For any two different sites 
i 7^ j, the creation and annihilation operators obey the 
usual bosonic relations [di, dj] = [dj, at] = [dj, a}-] =0. 

To understand the zero-temperature phase diagram of 
hardcore bosons in a superlattice potential, let us first 
analyze the atomic (t = 0) limit. In this limit, there is no 
kinetic (hopping) term, and the boson number operators 
hi commute with the Hamiltonian, so every lattice site 



is occupied by a fixed number of bosons. The average 
boson occupancy is determined so as to minimize the 
ground-state (free) energy. In particular, for A = 0, the 
model is translationally invariant, and the ground-state 
boson occupancy is the same for each of the lattice sites: 
for fi < the minimal energy configuration is simply the 
particle vacuum (VP), i.e., the completely-empty lattice, 
and for fi > the minimal energy configuration is simply 
the hole vacuum (VH), i.e., the completely- filled lattice. 
The ground-state energy of these phases is degenerate at 
[i = 0. When A ^ 0, the ground state has an additional 
half-filled insulating phase characterized by crystalline 
order in the form of staggered boson densities, i.e., (hi) = 
1 for the even (or odd, depending on the sign of fJ,/A) 
sublattice and (hi) =0 for the odd (or even) one. We call 
this alternating density pattern the MI phase, although it 
is sometimes referred to as a charge density wave.— The 
MI phase resides in the region \^/A\ < 1, sandwiched 
between the particle vacuum and the hole vacuum. 

Having discussed the t = limit, we are now ready 
to analyze the competition between the kinetic and the 
potential energy terms of the Hamiltonian when t ^ 0. 
In one dimension, the phase diagram of the model is al- 
ready known. As noted in the Introduction, the model in 
this case has an analytic solution! 20 ' 22 This is due to the 
Jordan- Wigner transformation which enables the map- 
ping of the hardcore boson Hamiltonian to that of non- 
interacting spinless fermions. The dispersion relation in 
this case is given by 

e(fc) = -ju ± yjAt 2 cos 2 (fca) + A 2 , (2) 

where a is the lattice constant. The phase diagram con- 
sists of three insulating incompressible regions (these are 
extensions of the t = ones), as shown in Fig. [Ha). Two 
are the 'trivial' insulators: the VP phase which is ob- 
tained for large and negative values of /i, and the VH 
phase which is obtained for large and positive values of 
U. These two phases are also present in the absence of 
the alternating potential, and are particle-hole 'mirror 
images' of each other. They are se parated from the SF 
phase along the curves /i/A = ±\/l + (2t/^4) 2 [see Fig. 
H(a)]. As evident from the expression for the dispersion 
relations given above, the superlattice (i.e., the onsite 
checkerboard potential) creates a gap of A = 2A in the 
energy spectrum, leading to a MI phase at half filling. 
This is the 'slab' enclosed by [if A = 1 from above and 
[i I A = — I from below, in the center of the figure. 

In dimensions higher than one [Fig.[]Tb)], the expected 
phase diagram of the hardcore Bosc-Hubbard model is 
qualitatively similar to the one-dimensional case with one 
notable exception. Here, the MI region does not extend 
to infinity, but instead is a finite lobe, connecting the two 
SF regimes together. 

The phase diagram of the hardcore Bosc-Hubbard 
model has one additional property resulting from it be- 
ing invariant under the transformation d 2 ; — > d| , . (where 
f denotes a shift of one lattice step in any of the possi- 
ble directions). This symmetry operation, which can be 
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FIG. 1: Phase diagram of the hardcore Bose-Hubbard model 
in the presence of a period- two superlattice, Eq. (fl]). In one 
dimension (left panel), the phase diagram contains three in- 
compressible insulating phases, indicated by VH - the hole 
vacuum, i.e., a completely-filled lattice, VP - the particle vac- 
uum, i.e., the completely-empty lattice, and MI - the Mott 
insulator, in which case the average density is 1/2 and the 
local densities on the even and odd sublattices are different. 
Outside of these insulating regions, the system is superfluid 
(SF). In higher dimensions (right panel), the phase diagram 
is similar, with one exception: while in one dimension the MI 
phase extends to infinity, in higher dimensions the MI phase 
takes the form of a Mott lobe. 



III. VACUUM OF PARTICLES AND HOLES 
PHASE BOUNDARIES 

As it turns out, the phase boundary separating the SF 
phase from the insulating VH phase (henceforth, the SF- 
VH boundary) can be easily obtained analytically for any 
given dimension. To sec this, we will use the fact that our 
Hamiltonian commutes with the total-number-of-bosons 
operator N = n,. In spin language, this simply means 
that for any given set of parameters /i, A and t, the 
ground-state wave function will be a linear combination 
of product states each having the same number of spin- 
downs. In the VH phase, this number is zero, as the wave 
function is simply 



|VH) = |ttt • • ■ ttt) , 



(5) 



with energy £ v h = —fiN. In the infinitcsimally thin layer 
outside the VH phase, the state of the system (which we 
shall refer to as the VH 'defect' state) is characterized by 
exactly one spin-down. That is, the wave function has 
the form: 



|VH dcf }= ]Tc4r |VH) 



(6) 



immediately read off from the Hamiltonian, corresponds 
to a particle-hole exchange combined with swapping the 
odd and even sublattices. It leads to a \i — > —fj, symme- 
try in the phase diagram. We shall make use of this fact 
when we obtain the phase diagram in later sections. The 
special case of /i = has been studied in Ref. Qjl both in 
two and three dimensions. 

Before moving on, we recall that the model at hand can 
also be viewed as the XY model of a spin- 1/2 system. 25 ' 26 
This is due to the mapping between bosonic operators 
and SU(2) generators: 



at o sf, 

ctj -f-j- Sj , 
a\cn o Sf + 1/2. 



(3) 



The symmetry of our model further tells us that all the 
coefficients Cj whose index V corresponds to a site on the 
even (odd) sublatticc arc all the same, namely: 



(7) 



where normalization requires 7V/2(|c c 



.I 2 ) 



= 1, 

and e.s. (o.s.) stands for the even (odd) sublattice. In 
order to determine the exact value of the weights c c s . 
and c .s., we first act with the Hamiltonian on this state. 
This eigenvalue problem then reduces to the following 
coupled equations: 

-2dt Co. s . + [m(1 -N) + A]ce. s . = e c e . s . (8a) 
-2dt Ce. s . + - N) - A}c . s . = e c Q . s . , (8b) 

where e is the energy of the state. Solving for e, the 
solution with minimal energy turns out to be 



With this mapping, the hardcore bosons Hamiltonian, 
Eq. ([1]) , becomes that of the XY antiferromagnet with an 
alternating magnetic field applied along the z direction: 



<'j) 



\<r(i) 



S ; 



(4) 



This alternative representation will become handy in the 
next sections. 



£dcf 



-fiN + \/A 2 + (2dty 



(9) 



The SF-VH boundary is the curve along which the VH 
state, Eq. ([5]), is no longer energetically favorable. This 
happens when its energy becomes equal to the energy of 
the defect state, Eq. ©. Matching the two, we obtain 
the SF-VH phase boundary: 



(10) 



where x — 2dt/A. 

A few remarks are now in order. As already noted in 
the previous section, the phase diagram of the hardcore 
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Bosc-Hubbard model is symmetric under the transfor- 
mation fj, — > — /Z. This tells us that the SF-VP phase 
boundary [the lowest branch in Fig. (Hb)], is given by 
fx/A = — y/l + x 2 . This result can also be obtained 
by repeating the above exercise with the substitution 
| t) «->• | I)- We also note that Eq. flTU|) agrees with 
the corresponding expression of the one-dimensional case 
obtained formerly (see Sec. IIIip . 

Another, simpler argument leading to the same solu- 
tion stems from the fact that the boundary between the 
SF and the VP (VH) phase is determined by the addi- 
tion of a single particle (hole) to the completely-empty 
(-filled) lattice. It can then be argued that whether one is 
dealing with hardcore bosons or noninteracting spinless 
fermions makes no difference in this the parti- 

cle statistics plays no role. This further means that one 
needs only to diagonalize the single-particle Hamiltonian 
and find the energy difference between the completely- 
empty (-filled) lattice and the state with one particle 
(hole). These will provide the chemical potential at the 
boundary between the SF and the VP (VH) phase. The 
single-particle spectrum in a d-dimensional superlattice 
with period two has the form: 

e{k) = -(jl± yjAdH 2 cos 2 (ka) + A 2 , (11) 
from which Eq. (|I0[) follows trivially. 



IV. NUMERICAL RESULTS 

Unlike the SF-VH and SF-VP phase boundaries, the 
SF-MI boundary, cannot be determined with the tools in- 
troduced in the previous section. One reason for that is 
that the exact many-body wave function of the MI state 
is not known. Therefore, in this section we explore the 
SF-MI phase boundary numerically by performing sim- 
ulations based on the stochastic series expansion (SSE) 
algorithm ! 23 ! 24 Our main objective here is to find the 
critical points of the SF-insulator transitions in the (x- 
A parameter space (without loss of generality we fix the 
hopping parameter at t — 1 and consider only fi > and 
A > 0). Critical points on the SF-MI boundary were typ- 
ically obtained by first fixing the value of the parameter 
A, and then performing the simulations for a range of val- 
ues of fi and different system sizes. This procedure was 
then repeated for different values of A. In some cases, 
mainly near the tip of the Mott lobe, we repeated the 
above procedure by fixing the value of fi and performing 
simulations for a range of A values and different system 
sizes. This was done mainly to further verify the accu- 
racy of the results, as the tip of the lobe is a multicritical 
point and therefore requires more care. 

Repeating the simulations with different system sizes, 
enables us to extrapolate the thermodynamic limit by 
correcting finite-size effects using scaling arguments in 
the vicinity of the phase transition: around the critical 
point, most physical quantities (which we denote here by 



X) scale according to the general rule: 

Xtf/v = F{\h-li c \L 1 I v ), (12) 

where F is a universal scaling function, fi — fi c is the 
shifted control parameter (/z being the control parame- 
ter, and fi c its critical value), v is the correlation length 
critical exponent and £ is the critical exponent belong- 
ing to the observable X . The values of these exponents 
are determined by the universality class the transition 
belongs to. In a previous workfi2. we studied the SF- 
MI transition at fixed (half-filled) density. This type of 
transition belongs to the (d+1) XY universality class, 
similarly to the SF to MI transition of the Bose-Hubbard 
model at fixed integer density^ Here, we compute the 
phase boundary between the SF and the (half-filled) MI 
phase while changing the density, so the transition be- 
longs to the mean-field universality class for which the 
correlation length and dynamical critical exponents are 
v = 1/2 and z = 2 (again, exactly as the corresponding 
transition in the Bose-Hubbard model) £ 

Equation (fT2"j) above will help us find the critical point, 
as it tells us that (a) the quantity XL>I V should be inde- 
pendent of the size of the system at the phase transition, 
and (b) when plotting XLr-l v against |/x — fi c \L x l v the 
resulting curve should be independent of the system-size 
as well. The quantity we shall be using to that end is 
the supcrfluid density, which has the critical exponent 
£ = v(d + z — 2) (see Ref. [H for details) where d is the 
dimension. 

We note here that since we are interested in the zero- 
temperature properties of the system, simulations are 
performed with high inverse-temperature (3 = 1/T (in 
our units, kg = 1), where in most cases we will find it 
sufficient to have (3 > 2L in order to obtain virtually zero- 
temperature results. (The effects of increasing (3 beyond 
this value are indiscernible.) 

As already discussed, in one dimension, our model has 
an analytic solution^ This enabled us to compare our 
numerical method against exact analytic results, as a 
check on our computational approach. No discrepancies 
between the analytical solution and the numerical one 
were found (see also Ref. 1 1 9f ) . 

In dimensions higher than one, no analytic solution to 
the model exists, so accurate results are obtainable only 
numerically. In the two dimensional case, we have ap- 
plied the SSE algorithm to systems of sizes ranging from 
16 x 16 to 48 x 48, with inverse-temperature j3 = 64. 
Figure [2] is an example of how scaling of the supcrfluid 
density data for the various system sizes is performed in 
order to find the critical point corresponding to A = 1.05. 
Here, the scaled supcrfluid density is plotted against fi 
for the different system sizes (the statistical errors of the 
quantum Monte Carlo simulations are on the order of 
magnitude of the symbol sizes). All curves intersect at 
/x c ps 0.178, signifying the phase transition for A = 1.05. 
The inset shows the scaled superfluid density as a func- 
tion of the scaled control parameter, in which case all 
curves should be, and in fact are, on top of each other. 
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The resulting SF-MI phase boundary of our model in 
two-dimensions is marked by the full circles in Fig. [3] As 
noted earlier, the lower half of the phase diagram Fig. 
[T] (the fj, < half) is but a mirror image of the portion 
shown in Fig. [3j and thus is not presented there. The 
tip of the Mott lobe was found to be at x c ~ 2.Q2} 9 
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FIG. 2: (Color online) Scaled superfluid density as a function 
of the chemical potential /i for the various system sizes in 
the two-dimensional case (here, A — 1.05). All the curves 
intersect at fi ~ 0.178 indicating the value of the critical point. 
In the inset, the control parameter (the horizontal axis) is 
scaled as well, leading to the collapse of all data points into a 
single curve. 

In three dimensions, we have performed simulations 
with system sizes ranging from 6x6x6 to 16x16x16 
and an inverse temperature of j3 — 40. Figure U is an 
example of how scaling is carried out in three dimensions: 
the scaled superfluid density is plotted as a function of 
fi for the different system sizes and A = 2.28. The inset 
depicts the scaled superfluid density as a function of the 
scaled control parameter, exhibiting the collapse of all 
data points into a single curve, as in two dimensions. The 
resulting phase boundary in three-dimensions is shown in 
Fig. E] (full circles). The tip of the Mott lobe was found 
to be at x c « 1.44^2, 

V. MEAN-FIELD APPROACHES 

Having obtained the exact boundaries of the phase di- 
agram of the model, we now proceed to study several ap- 
proximation schemes, and examine the extent to which 
they provide an accurate description of the phase dia- 
gram of the model. We start this investigation with the 
Gutzwiller mean-field approach. 
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FIG. 3: (Color online) Phase diagram of the model in two di- 
mensions. The full circles are the analytical (VH boundary) 
and numerical (MI boundary) results. The solid line corre- 
sponds to the strong-coupling expansion (SCE) fit, whereas 
the dot-dashed and dashed lines are the mean-field (with and 
without spin-wave corrections) and cluster mean-field predic- 
tions, respectively. As the figure shows, the SF-VH bound- 
ary is predicted correctly by the mean-field approximation 
schemes. As for the SF-MI boundary, the predictions of the 
SCE fit provide the most accurate results. 
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FIG. 4: (Color online) Scaled superfluid density as a function 
of the chemical potential fi for the various system sizes in 
the three-dimensional case (here, A = 2.28). All the curves 
intersect at /i ~ 0.752 indicating the value of the critical point. 
In the inset, the control parameter (the horizontal axis) is 
scaled as well, leading to the collapse of all data points into a 
single curve. 



6 



< 

~5 



0.5 




FIG. 5: (Color online) Phase diagram of the model in three 
dimensions. The full circles are the analytical (VH boundary) 
and numerical (MI boundary) results. The solid line corre- 
sponds to the strong-coupling expansion (SCE) fit, whereas 
the dot-dashed and dashed lines are the mean-field (with and 
without spin-wave corrections) and cluster mean-field predic- 
tions, respectively. As the figure shows, the SF-VH bound- 
ary is predicted correctly by the mean-field approximation 
schemes. As for the SF-MI boundary, the predictions of the 
SCE fit provide the most accurate results. 



A. Gutzwiller mean-field 

Along the lines of Ref. HH, we start our mean-field cal- 
culation with the following product state as our ansatz: 



10} 



MF 



Ijsinf ||)+cos^e 



■<Pi I 



t) 



(13) 



The angles (9j, ipj) here, specify the orientation of the j- 
th spin. Naturally, we expect the wave functions of each 
of the odd (even) sublattice sites to be identical. This is 
due to the checkerboard symmetry of the model. 

As we are using the grand-canonical scheme, the ori- 
entations of the spins will be determined by minimizing 
the grand-canonical potential (per site) 
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with respect to these angles. For the azimuthal an- 
gles, this simply implies a constant (yet arbitrary) value 
(fj = while for the polar angles, the minimizers are 



cos 9± = Min 
cos 6*2 = Min 
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where /ii j = (p ± A)/(2dt). We note that while in Ref. 
IT9I the focus was on the special p = case, here we place 
no limitations on p. 

At this point we can calculate the following quantities. 
First, the average density of particles is: 
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Next, the free energy becomes 

^MF = MF(0|-ff|0)MF = 

(p + A) cos 6»i 



MF \ 
1 

4 



dt . n . n 
— — sin 9\ sin 62 — 



{p~A)cosd 2 , (17) 



and the density of bosons in the zero-momentum mode 
po is calculated as: 
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(18) 



The superfluid density too is obtained in a straight- 
forward manner. In the mean-field approximation it has 
the simple form p s = — (2d)~ 1 dfl/dt.— 

The phase boundaries are simply the curves along 
which the superfluid density and the zero-momentum 
fraction drop to zero. These turn out to be: 



(19) 



where the '+' branch belongs to the SF-VH transi- 
tion and the '-' branch belongs to the SF-MI transition 
(again, x = 2dt/A). The phase diagram of the model 
as predicted by the Gutzwiller mean-field approach is 
sketched in Fig. [SJ which shows the average density of 
bosons as a function of x and pj 'A. 

An alternative way of deriving the mean-field phase 
boundaries is through the decoupling approximation.^! 
In this approach, one approximates the hopping term as 



a- a,4 
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(20) 



and introduces the condensate order parameter ipi = 
sfrii = (aj) = (a,i) (analogous to the Bogoliubov ap- 
proach). Since the condensate order parameter is the 
same for all lattice sites belonging to the same sublat- 
tice, i.e., 



Ipc.s. + Ipo.t 



+ (-1) 



<r(i) 



(21) 



for some real unknown parameters ipe.s. and ipo.s. (due to 
the checkerboard symmetry of the model) , it is sufficient 
to solve only for the effective two-site Hamiltonian 
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FIG. 6: (Color online) Average density of bosons as a function 
of x — 2dt/A and n/A in the mean-field approximation. The 
three insulating phases VP (empty lattice, zero density), MI 
(half-filled lattice) and VH (completely-filled lattice) are seen 
very clearly in the figure. Outside of these insulating regions 
is the SF phase. 



where i € e.s. and j € o.s.. Performing a second-order 
perturbation theory in the first two terms of this effec- 
tive Hamiltonian around the VH and MI phases produces 
the ground state energies as a function of ipe.s. and V'o.s.- 
Notice that higher orders are not needed for our pur- 
poses, since the second order theory is sufficient to derive 
the energy functional of the system up to second order 
in the order parameters ipo.s. and ip c .s.- Following the 
usual Landau procedure for second-order phase transi- 
tions, minimizing the ground state energies as a function 
of the superfluid order parameters, we eventually arrive 
at Eq. (PI). 

The dash-dotted lines in Figs. [3] and [5] show the phase 
diagram as predicted by the Gutzwiller mean-field ap- 
proach, compared against the numerical results. Inter- 
estingly, the mean-field ansatz yields the correct solution 
for the SF-VH transition (upper branch). On the other 
hand, for the SF-MI boundary, mean-field results differ 
considerably from the numerical data: while away from 
the tip of the Mott lobe the method is very accurate, 
as one approaches the tip itself, errors climb up to their 
maximal values of ~ 100% in two-dimensions and ~ 50% 
in three dimensions at the tip of the MI lobe. The very 
large errors here reflect the fact that the mean-field ap- 
proach is simply not fit to describe the phase transition 
in this region. 

Before moving on, we remark here that addition of 
spin-wave corrections to the mean-field solution does not 
modify the mean-field critical points of the model^ so 
the phase boundary is not altered by spin-wave correc- 
tions. While deep in the SF phase spin-wave corrections 
yield major improvements over the mean-field results for 
many of our obscrvablcs of interest, in terms of phase 



boundaries the spin-wave corrections do not contribute. 
As one approaches the phase transition itself, the spin- 
wave corrections lose their accuracy, eventually leaving 
the phase boundaries at their mean-field values 

B. Cluster mean-field 

Aiming to improve the results obtained in the previous 
section, we now describe a 'cluster' mean-field approach, 
which makes use of the checkerboard symmetry of the 
model. This approximation scheme was introduced in 
Ref. Qj^ where it was applied to the special case of /i = 
0. Within this approach, one starts with a variational 
ansatz which, as before, is a product state. However, 
this time one docs not choose a product of single-site 
wave functions. The new ansatz is a product of wave 
functions each describing the state of a 'block' of 2 d sites, 
such that with this block as the basic cell, the model turns 
homogeneous. In two dimensions, for example, a block 
consists of 2 x 2 square cells each of which is described 
by the general wave function 

|0)cmf= J] f Yl c w\ijkl)) , (23) 

blocks \i,j,k,le{lA} J 

where the generalization to three dimensions, in which 
case the basic block is a 2 x 2 x 2 cubic cell, is straight- 
forward (note that the coefficients for each of the blocks 
will be the same due to the symmetry of the wave 
function)^ As before, we minimize the free energy 
^CMF = CMF(0|-ff |0)cmf with respect to the coefficients 
Cijki of the wave function (this time we do so numeri- 
cally) . Obtaining the various observables in terms of the 
wave function given in Eq. (|23[) is straightforward, and 
was performed in much the same way as the usual mean- 
field approach discussed in Sec. IV A[ 

The phase boundaries, as predicted by the cluster 
mean-field approximation, are given by the dashed lines 
in Figs. |3] and [5] for two and three dimensions, respec- 
tively. As the figures indicate, the SF-VH boundary is 
predicted correctly. This is no surprise as the Gutzwiller 
mean-field, over which the current method is an improve- 
ment, is already exact for that boundary. As for the SF- 
MI boundary, the cluster mean-field method is far better 
than the Gutzwiller mean-field method. As in the previ- 
ous mean-field case, the results are more accurate away 
from the tip of the Mott lobe but reach ss 60% error in 
two dimensions w 24% error in three dimensions, as the 
tip is approached. 

Having shown that the mean-field-type theories pre- 
sented here are not very accurate in describing the SF- 
MI phase boundary, in particular close to the tip of the 
lobe, we turn to develop a strong-coupling perturbation 
theory in the hopping t. This approach, combined with a 
scaling analysis, will allow us to predict the critical point 
and the shape of the insulating lobe in a more accurate 
manner. 
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VI. STRONG-COUPLING EXPANSION (SCE) 

Strong-coupling expansion (SCE) techniques were pre- 
viously used to discuss the phase diagram of the Bose- 
Hubbard model^^^ and of the extended Bosc-Hubbard 
models and its results showed an excellent agreement 
with quantum Monte Carlo simulation o n i 12 in the for- 
mer case. Motivated by the success of this technique 
with Bosc-Hubbard type models, here we generalize this 
technique to the hardcore Bose-Hubbard model on a su- 
pcrlattice. 

To determine the phase boundary separating the in- 
compressible MI phase from the compressible SF phase 
within the SCE method, one needs the energy of the MI 
phase and its 'defect' states - those states which have one 
flipped spin (equivalently, one excited particle) about the 
ground-state - as a function of the parameter t. At the 
point where the energy of the incompressible state be- 
comes equal to its defect state, the system becomes com- 
pressible, assuming that the compressibility approaches 
zero continuously at the phase boundary. Note that these 
arguments are very similar to those presented in Sec. IIIII 
where exact results were obtained for the SF-vacuum in- 
sulators boundaries. Here however, the state of the sys- 
tem inside the MI phase is not known except for the 
special case t = 0, where: 

|mi<°>) = inn •••nn), (24) 

where all the spin-ups (spin-downs) belong to the even 
(odd) sublattice. 

The energy of the MI phase is calculated via a many- 
body version of the nondegenerate Rayleigh-Schrddingcr 
perturbation theory up to fourth order in t. We note 
that all odd-order terms in t vanish for the d-dimensional 
hypercubic lattices considered in this manuscript. This 
is because this state cannot be connected to itself by 
only one hopping, but rather requires two hoppings to 
be connected. 

Calculation of the wave functions and energies for 
the defect states is more involved as it requires the use 
of the many-body version of the degenerate Rayleigh- 
Schrodinger perturbation theory. The reason for that 
lies in the fact that when exactly one extra particle is 
added to the MI phase, it could go to any of the N/2 
lattice sites that belong to the odd sublattice, since all 
of those states share the same energy when t = (recall 
that N is the number of lattice sites). Therefore, the 
initial degeneracy of the MI defect state is of order N/2. 
This degeneracy is lifted at second order in t, since all of 
the defect states occupy one of the sublattices, and they 
cannot be connected by one hopping, but rather require 
two hoppings to be connected. The wave function (to 
zcroth order in t) of the particle-defect state turns out to 
be 

|MO = J2 M^|MI(°>), (25) 



where /, is the eigenvector of the matrix T^i = 
52jee s tijtji' with the highest eigenvalue, such that 

E.'eo.sT"'/!' = 4d 2 i 2 / l . Here, Uj = t for (ij) and 
zero otherwise. The normalization condition requires 
that J2ieos \fi\ 2 = 1- The eigenvector with the highest 
eigenvalue corresponds to the lowest energy state, i.e., to 
the ground state. We calculate the energy of the |MI*?l) 
phase via degenerate perturbation theory up to fourth 
order in t. Here too all odd-order terms in t vanish. 

A lengthy but straightforward calculation leads to the 
following expression for the SF-MI boundary (for further 
details regarding the calculation, we refer the reader to 
a similar calculation given in Ref. [29t ) 



£ , a (d-l)(d-3) 4 

A 2d X 8d 2 X 



where x = 2dt/A. This expression is exact for all d- 
dimensional hypercubic lattices up to the given order. In 
one dimension, it agrees with the analytical solution 2 ^ of 
the model given by \i/A = 1 (see Sec. UTJ) . In the d — >• oo 
limit, where the exact result is given by the mean- field 
expression, i.e., fi/A = \/l — x 2 , Eq. (|26|) is the correct 
power-series expansion about x = 0. 

In the two- and three-dimensional cases, fourth-order 
SCE is not very accurate near the tip of the MI lobe, as 
the variable x is not very small there. Therefore, an ex- 
trapolation technique is desirable in order to determine 
the phase boundary more accurately. Such an extrap- 
olation is possible for the MI phase, since it is already 
known for d > 1 that the critical point at the tip of the 
MI lobe has the scaling behavior of a (d+1) XY model. 
Therefore, we propose the following ansatz for the MI 
lobe which includes the known power-law critical behav- 
ior of the tip of the lobe: 

— = Qo (1 + oiix + a^x + a$x + ot^x ) 

x {x c -x) zv , (27) 

where x c — 2dt/A c is the critical point which determines 
the location of the MI lobe tip, and zv is the critical 
exponent for the XY model which determines the 

shape of the MI lobe near x c . The parameters oii are 
determined by matching Eq. (|26[) with Eq. (|2"T|) . after 
the latter is expanded out to fourth order in t. This 
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procedure leads to: 



a 

Oil 



x\ 
z v 



1 

ZV 



+ e 2 , 



zv{zv + 1) 

zv(zv + l)(zv + 2) 
h — e2 ' 

zv(zv + l){zv + 2)(zv + 3) 



+ 1) 
2xJ 



-e 2 + e 4 , 



(28a) 
(28b) 
(28c) 
(28d) 

(28e) 



where e 2 = -(d-l)/(2d) and e 4 = -(d- l)(d-3)/(8d 2 ) 
are the coefficients of the second and fourth order terms 
in our SCE. 

In our extrapolations, we set zv m 0.672 for d = 2 
and zv = 1/2 for d > 2. This leaves only x c to be fixed; 
something which is accomplished by a straightforward x 2 
curve-fitting to the numerical data obtained in Sec. IIVI 
The results are shown by the solid lines in Figs. [3] and [5] 
for two and three dimensions, respectively. As one can 
immediately see, the SCE results are very accurate and 
provide an analytic expression for the phase boundaries. 

Alternatively, we can estimate x c using the above ap- 
proach without fitting it to the numerical data. We do so 
by finding the value of x c for which the fifth-order term 
in x of Eq. (|27p vanishes. This gives x c « 1.53 for d = 3 
(« 6.7% error), and x c « 1.076 for d -> oo (w 7.6% 
error). 

Before moving on to the next section, we note here that 
a similar application of the SCE for the SF-VH phase 
boundary, where 



£ MrivH) 



(29) 



? Go.s. 



is the wave function (to zeroth order in t) of the hole- 
defect state, leads to 



A = 1+ r 



i 



-a; 4 + 0(a; 6 ), 



(30) 



in agreement with the exact result derived in Sec. lIIIl i.e., 
fj,/A = yl + x 2 , up to the given order. In addition, we 
perform a SCE in A, and find that the large x behavior 
of the phase boundary is given by fi/A = x + 0(l/x), 
which is also in agreement with the exact result . 



VII. MOMENTUM DISTRIBUTION 

Having discussed the phase diagram of the hardcore 
Bosc-Hubbard model with a supcrlatticc in the previ- 
ous sections, next we analyze the momentum distribu- 
tion n(k) of these bosons. This quantity can be directly 



probed in experiments with ultracold atomic gases via 
an absorption imaging during a short time-of- flight 
Since it is trivial to show that nvii(k) = 1 in the VH 
phase, we shall concentrate only on the momentum dis- 
tribution of the bosons in the MI phase, riMi(k), where we 
will compare our numerical quantum Monte Carlo results 
with those of two analytical approaches: the random- 
phase approximation (RPA) and the SCE method intro- 
duced in the previous section. 

The RPA is a well-defined linear operation in which 
thermal averages of products of operators are replaced by 
the product of their thermal averages^ Since the fluctu- 
ations are not fully taken into account in this method, it 
becomes exact only for infinite-dimensional bosonic sys- 
tems, recovering the mean-field theory. This method has 
been recently applied to the onsite ; 32 ' 33 and extended 3 ^ 
Bose-Hubbard models, and its results showed good qual- 
itative agreement with the experiments in the former 
c&se££2- Here we apply this method to our model (for fur- 
ther details regarding the calculation, we refer the reader 
to a similar calculation given in Ref. I34I). and obtain 



"MI-RPA(£k) — - 



lA- 




A + 





(31) 



where = ~ 2t 2<=i cos(/cjo) is the energy dispersion of 
noninteracting particles. Since the RPA phase boundary 
is exactly the same as the mean-field one, and it gives a 
critical value for x = 2dt/A that is much smaller than 
the true critical value in finite-dimensions, we compare 
our results with a rescaled x value such that 



scaled / \ 

"MI-RPA^kJ — 



1 Ax 



£k 



Ax 



£k 



(32) 



where x c = 2dt/A c is the true critical point which deter- 
mines the location of the MI lobe tip. We call this the 
scaled RPA momentum distribution following Ref. H^. 

To extend the RPA result to finite dimensions, we also 
calculate nMi(ek) as a power series expansion in the hop- 
ping t via the strong-coupling perturbation theory. To 
second-order in t, we obtain (for further details regard- 
ing the calculation, we a gain refer the reader to a similar 
calculation given in Ref. KM ) 



, v 1 £k 

"-MI-SCE(£kJ - - - 2^ 



4 - 2dt 2 
AA 2 



+ 0(t 3 ), (33) 



which is exact up to the given order for any dimension d. 
In the d — > oo limit (while dt is kept fixed), we checked 
that Eq. ([3"3"]) agrees with the RPA solution (which is 
exact in that limit) given in Eq. (|31[) . when the latter 
is expanded out to second order in t. This provides an 
independent check of the algebra. 

The second-order SCE is not very accurate near the 
tip of the MI lobe, as t/A is not small there. To extend 
its region of validity, we therefore propose the following 
ansatz, 



. 1 A-e k + (4X-2)dt 2 /A 
nMl(£k) = 2\/ AT^+4Xa¥jA {M) 
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FIG. 7: (Color online) Momentum distribution function n(ek) for two 48 x 48 systems: (a) A = 2.42 (a; ~ 1.653) and (c) 
A = 4.22 (x w 0.948), and two 14 x 14 x 14 systems: (b) A = 4.8 [x = 1.25) and (d) A = 6 (x = 1). The full circles are the 
numerical SSE results. The scaled RPA is the dot-dashed line, whereas the dashed and solid lines are the SCE and scaled SCE, 
respectively. The figures show that the scaled SCE results are much better than any of the other two approximation methods, 
and that the scaled SCE results fit better, as A becomes larger (t = 1 in all four systems) - suggesting we are deeper inside 
the MI phase. 



for any dimension d, where A = d(x c — l)/x^ depends 
on d. This expression reduces to Eq. (|3~T|) in the d — > oo 
limit, and it has the correct power-series expansion about 
x = up to second-order in t, i.e., Eq. We call this 

the scaled SCE momentum distribution. 

In Fig. we show several comparisons (two in two 
dimensions and two in three dimensions) between the 
momentum distribution function obtained with the quan- 
tum Monte Carlo and the three approximations obtained 
above, namely, the scaled RPA, the SCE, and the scaled 
SCE. As the figures indicate, the scaled SCE is a far bet- 
ter fit than the other two methods, and more so for larger 
values of A, that is, deeper inside the MI phase where the 
SCE becomes more and more accurate. 



VIII. CONCLUSIONS 

We have obtained the complete phase diagram of 
the hardcore Bosc-Hubbard model with a period-two 
supcrlatticc in two and three dimensions. First wc 
have calculated the boundaries between the superfluid 



phase and the 'trivial' insulators (the completely-empty 
and completely-filled lattices) analytically. Then, us- 
ing quantum Monte Carlo simulations followed by a 
finite-size scaling, we have determined the phase bound- 
ary between the superfluid phase and the (half-filled) 
Mott insulator. We have also compared our numer- 
ical results against three approximation schemes: the 
usual Gutzwiller mean-field approach, a cluster mean- 
field approach, and the strong-coupling expansion (SCE) 
method. 

For the transition between the superfluid phase and 
the 'trivial' completely-empty and completely-filled lat- 
tice insulators, we have found that the mean-field ap- 
proaches yield the exact results in any dimension. As for 
the supcrfluid-Mott insulator boundary, the Gutzwiller 
approach was shown to work very poorly (up to 100% 
error in two dimensions and ~ 50% error in three di- 
mensions). This is a clear indication of the fact that 
this mean-field approach is not suitable for describing 
the superfluid-Mott insulator transition in the vicinity of 
the tip of the lobe. A cluster mean-field approximation 
scheme, which is based on the underlying checkerboard 
symmetry of the problem, was proven to be a big im- 
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provement over the previous method (reducing the error 
to one half of the one generated by the usual Gutzwiller 
ansatz), albeit still far from being accurate as one ap- 
proaches the tip of the Mott lobe. The fourth-order SCE 
turned out to be the best method among the three in de- 
scribing the superfluid-Mott insulator phase boundary, 
as the one-parametric fit of the SCE yielded very accu- 
rate results, also near the tip of the Mott lobe where the 
other methods failed. It also provided an analytic expres- 
sion for that boundary, which could be used as a guide 
in future experimental realizations of this model. 

Finally we have examined the extent to which several 
approximation schemes, such as the random phase ap- 
proximation and the strong-coupling expansion, give an 
accurate description of the momentum distribution of the 
bosons inside the insulating phases. We have shown that 



a scaled SCE provides an accurate analytic expression 
for the momentum distribution of the bosons inside the 
Mott-insulating phase both in two and three dimensions, 
which could again be used as a guide in future experi- 
mental realizations of this model. 
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